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Abstract 

Based on a high temperature expansion, we compute the two-point correla- 
tion function and the critical line of an Ising lattice gas driven into a non- 
equilibrium steady state by a uniform bias E. The lowest nontrivial order 
already reproduces the key features, i.e., the discontinuity singularity of the 
structure factor and the (qualitative) ^'-dependence of the critical line. Our 
approach is easily generalized to other non-equilibrium lattice models and 
provides a simple analytic tool for the study of the high temperature phase 
and its boundaries. 

I. INTRODUCTION 

The study of non-equilibrium steady states (NESS) has attracted vivid interest over the 
past decade. On the one hand, such studies are application-driven, since NESS determine 
the physics of a wide range of important problems, including, e.g., granular and traffic flow, 
surface growth, electromigration, and transport phenomena in biological systems. On the 
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other hand, there is a fundamental interest in creating a theoretical framework for NESS on a 
par with Gibbs ensemble theory for equilibrium systems. Driven diffusive lattice gases, first 
introduced by Katz, et.al. and recently reviewed in 0, provide simple testing grounds 
for the properties of a particular class of NESS. Based mostly on Ising lattice gases with 
Kawasaki (spin exchange) dynamics, these models are represented by microscopic Master 
equations which violate the usual detailed balance condition by maintaining a net probability 
current between configurations. A particularly simple member of this class is the uniformly 
driven lattice gas, in which spin exchanges along a specific lattice axis are biased by a uniform 
force, remaining energetically controlled only in the transverse subspace. Identifying the 
particles with (shielded) charges and the force with an electric field, such models provide 
a starting point for the description of fast ionic conductors or charged droplets in a 
microemulsion [Q. 

One of the more remarkable features of this model, and many of its relatives, is the 
presence of long-range correlations at all temperatures Though this behavior can be 

understood, within the context of a phenomenological (mean) field theory in terms of the 
violation of the traditional fluctuation-dissipation theorem it is clearly important to have 
an exact microscopic version. Using a venerable tool, the high temperature series expansion, 
Zhang, et.al. P] (ZWLV in the following) investigated G, the two-point correlation function 
for a system with infinite E. From the Master equation and the associated hierarchy, they 
argued that three-point correlations are negligible and arrived at a closed set of equations for 
G alone. To zeroth order in (3, both the equation and the solution are trivial. To first order, 
the short- distance behavior of the solution was obtained numerically, by truncating the 
equations at distances larger than a cutoff value. The results exhibit significant anisotropy 
and agree quite well with Monte Carlo data at high temperatures, i.e., T > 6, in units 
of Tc(0), the critical temperature of the E = system. Analytically, the equations were 
approximated by a Poisson problem with quadrupole symmetry. Its solution captures the 
behavior of G for large distances, displaying the power law tail, with the angle-dependent 
amplitude of the quadrupole potential, in agreement with simulations. We should emphasize 
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that the high temperature series has a strong foundation: for (3J = 0, the steady state 
distribution P* oc 1 is exactly known for all ||^, so that we are expanding about a well- 
defined zeroth order state. In this paper, we generalize the analysis of ZWLV to include 
finite driving fields and solve the resulting set of equations for G exactly, by computing its 
Fourier transform, the structure factor S. The latter displays the characteristic discontinuity 
singularity at the origin, which translates into the r"'' decay in real (rf-dimensional) 
space. Thus, recourse to numerical methods is not necessary. 

Another key feature of driven lattice gases is the existence of a line of continuous tran- 
sitions, falling outside the Ising universality class [0. In particular, it is remarkable that 
the critical temperature, Tc{E), increases with E, even though the drive reduces the effec- 
tive nearest-neighbor coupling. Apparently, the strong long-range part of the correlations 



dominates the driven system and permits the onset of order at a higher 1TT[]. To support 
this argument, it is desirable to have a simple analytic method to estimate both correlations 
and the critical temperature. The high temperature expansion, and specifically our explicit 
form for the structure factor, serve that purpose. Matching the series expansion of S to the 
expected critical singularity, we obtain estimates for the critical line, and hence the phase 
diagram to this (lowest) order in /3. 

Clearly, quantitative accuracy cannot be expected from a first-order calculation. Never- 
theless, it presents one of the simplest, currently available, tools for the qualitative prediction 
of nonuniversal quantities in driven lattice gases. Real-space renormalization group tech- 
niques for conserved non-equilibrium spin systems are sorely lacking, and dynamic mean- 



field theories [|T2|, while quantitatively more satisfying, are much more labor-intensive. Our 



method is easily extended to general rate functions, anisotropic interactions and higher 
dimensions. These results will be published elsewhere |I3| . 

The paper is organized as follows: we first summarize the model definition and the 
method of ZWLV, resulting in a closed set of equations for the two-point function G, correct 
to first order in (3J but for arbitrary drive strength E > 12 J . These equations are then 
solved exactly in Fourier space. Key consequences, such as the discontinuity singularity of 



the structure factor and the critical hne Tc{E), are discussed. We conclude with a brief 
comment on results in higher dimensions. 



II. THE MODEL AND THE EQUATIONS FOR ITS TWO-POINT FUNCTION 

We first summarize the microscopies of the model. On each site r of an (infinite) square 
lattice in spatial dimension d = 2, coupled to a heat bath at inverse temperature f3, we 
define an Ising spin variable 5*^ which distinguishes occupied {S = +1) from empty (5* = —1) 
sites. Interacting via the usual Ising Hamiltonian, Ti = — J^Sf^Sf^, nearest neighbors can 
exchange positions, subject to the local energetics and a uniform force E ("electric field") 
which biases exchanges along a specific lattice axis (aligned with the a;-coordinate) . To be 
specific, we choose Metropolis [1^ rates, min{l, exp —P[AH — aE]}, where a = 0, +1, —1 for 



jumps transverse, along and against the field. Thus, "infinite" E implies that jumps along 
the field always take place, provided they are allowed by the excluded volume constraint, 
while jumps against E are completely suppressed. These rates specify the master equation of 
the model, for the time-dependent configurational probability P. A hierarchy of equations, 
connecting different A^-point functions, follows as usual. We will be interested in steady- 
state averages only, so all time derivatives will be set to zero. Following ZWLV, we expand 
the rate functions in powers of (3 J; as an extension of their work, which focused on = cxo 
only, we include the effect of large, but finite (3E. Thus, in a sense, the case here corresponds 
to (3E ^ (3J . In practice, it is sufficient to choose E > 12 J, so that jumps along the drive 
occur with unit rate, while those against E are suppressed by a factor of exp{—j3E). In 
a strict high temperature expansion, organized in powers of f3, (3E would also appear as a 
small parameter, so that the range E < 12 J can also be explored |13[. 

Returning to our case, a set of equations for the two-point function, G{f, P) =< S^Sj^ >, 
is easily derived from the master equation with the expanded rates. At first order in (3 J, 
one finds that G couples only to three-point functions. Even though, in contrast to the 
equilibrium Ising model, the latter do not vanish here []T3[, they are numerically quite small 
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and are neglected. A closed set of equations for G emerges. By translational invariance, G 
depends only on {x,y), the separation between the two points f and f", and it is even in 
both variables. 

The resulting equations depend on two parameters, 

K = PJ and e = exp(-/?£;) . 

Up to and including only first order terms in we obtain a set of linear equations for G, 

= 2[G'(1, 1) + ^(1, -1) - 2^(1, 0)] + (1 + e)[G{2, 0) - ^(1, 0)] + 2i^(2 + e) 

= 2[G(0, 2) - G(0, 1)] + (1 + e)[G(l, 1) + G(-l, 1) - 2G(0, 1)] + 2K{1 + 2e) 

= 2[G(1, 2) + G(l, 0) - 2G(1, 1)] + (1 + e)[G(2, 1) + G(0, 1) - 2G(1, 1)] - 2K{1 + e) (1) 

= 2[G(2, 1) + G(2, -1) - 2G(2, 0)] + (1 + e)[G(3, 0) + G(l, 0) - 2G(2, 0)] - 2eK 

- 2[G(0, 3) + G(0, 1) - 2G(0, 2)] + (1 + e)[G(l, 2) + G(-l, 2) - 2G(0, 2)] - 2K 

and, for all other non-zero x, y: 

= 2[G{x, y + 1) + G{x, y-l)- 2G{x, y)] 

+(1 + e) [G{x + 1, y) + G(x - 1, y) - 2G{x, y)] (2) 

Here, we have written the equations in such a form that the left hand side would simply be 
dtG{x,y). The ZWLV equations can be retrieved by setting e = 0. The solution to these 
equations is most easily found in Fourier space. In other words, we consider the structure 
factor: 

oo 

S(k,p) = ^G(x,y) ex.p[-i(kx + py)] 

which is real, since G is even. Some details on Fourier transforms are given in Appendix A. 
In the absence of J, the solution is trivial: 

G(0,0) = 1 
G{x,y)^0 x,y^O, 
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so that 

S{k,p) = 1 = S 

reflecting, of course, that P* oc 1 in this case 0. So, the information about interactions is 
carried by S, defined by: 

s = s + s . 

The correlations appearing in being proportional to K, are really the transforms of S: 
G{x,y)= / S{k,p) expi{kx+py) , (x, 7^ (0, 0) (3) 



For convenience, we have introduced the notation 
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dk / dp 



(27r)2 

and will denote the anisotropic lattice Laplacian by 

A{k,p) = 2(1 + e)(l - cos A;) +4(1 - cosp) . 

Inserting (|^) into (|T], 0) and using the fact that 5* is real, we obtain 

J SA exp{ik) + (1 + e) J 5(1 - cos k) = +2(2 + e)K 
j SAex-p{ip) + 2 j S{1- cosp) = +2(1 + 2e)K 

J SA exp[i{k + p)] = -2(1 + e)K (4) 
j SAexp{2ik) = -2eK 
j SAexp{2ip) = -2K 

with similar equations for negative values of a; or 7/ and 

j SA ex'p[i{kx + py)] =0 for |a;| + |?/| > 2 . (5) 

Next, we seek to invoke the completeness relation (see the Appendix) in order to project 
out an equation for S. While the two additional terms on the left hand side of @) might 



appear to spoil this approach, we can just treat them, for the time being, as unknown 
e-dependent coefficients, 

h = j 5(1 - cos fc) 

S{1 — cosp) (6) 



and move them to the right hand side. Finally, for completeness, we need an additional 
equation for x = y = 0, namely: 

y"^A = [2(l + e)(l -cosfc) +4(1 -cosp)] = 2(1 + e)Ji +4/2 . (7) 

Now that we have equations for a//integer values of x, y, we can use ^ exp[i{kx+py)] = 
{2nY5{k)5{p). The result is 

S{k,p)=L{k,p)/A{k,p) (8) 

where L is the sum of terms on the right hand side, i.e., 

L{k,p) = 2/i(l + e)(l - cos A;) + 4/2(1 - cosp) 

+4/r[e(l -cosA;) + (1 -cosp)] [l + 2cosA; + 2cosp] . (9) 

Of course, (§) is still an implicit equation for the structure factor, due to the appearance 
of /i and I2 on the right. To find the explicit solution, we must determine the /'s. Since we 
have two unknowns, we need two linearly independent equations. The first of these follows 
from (0): Inserting our result for L{k,p), Eqn. (Pf) into the ffist equation of (H), we obtain 



Hk,p) 

A(A:,p)' 



= -/i + / ^,,;-^; (l - cos k) = M^jlj + KNi 



where 



Mi,i = -1 + 2(1 + e) y (1 - cosky/A{k,p) 
Mi,2 = 4 y"(l - cosA;)(l - cosp)/A(fc,p) 



and 



iVi = y (l-cosA;) [4(1 - cosp) +4e(l - cosA;)] [2 cos A; + 2 cosp + 1] /A{k,p) 
Due to a remarkable symmetry under {k,p) exchange [|13|, the second equation in (P) is 



not hnearly independent from the first. Instead, an additional equation is provided by the 
value of G at the origin, i.e., 1 = G{0, 0) = / S{k,p) which leads to 



= ls{k,p) = I = + KN, 



and 



with 

M2,i = 2(1 + e) cosk)/A{k,p) 
^2,2 = 4 j {l-cosp)/A{k,p) 

N2 = j [4(1 - cosp) +4e(l - cosA;)] [2 cos A; + 2 cosp + 1] /A{k,p) . 

Note that the singularity of l/A(A;,p) at the origin is cancelled by zeros in the numerators, 
so that all of these integrals are perfectly finite. The explicit solution now follows as 

= -KM'l^Nn . (10) 

Together with (^, this determines the full structure factor, displaying the expected propor- 
tionality S{k,p) oc K . 

For later reference, let us briefly consider the equilibrium analog to our results so far. 
To first non-trivial order in K, the two-point correlation of the Ising model is given by the 
well-known form [^ 

Ge9(0,0) = 1 

G,,{1,0)=G,,{0,1)=K 
Geq{x,y) = for all other x,y . 

This results in a structure factor Seq{k,p) = l-|-2_ft'(cos A;-|-cosj9) which is of course isotropic. 



III. DISCONTINUITY SINGULARITY AND CRITICALITY 



In this section, we discuss two of the most interesting consequences of the full solution, 

S = S + ^ + OiK'') . (11) 



First, we focus on the celebrated discontinuity of the structure factor near the origin [0,|1|. 
This anomaly is a direct consequence of FDT violation and is therefore expected to increase 
in magnitude with the strength of the drive. For small k,p we find 

A = (1 + e)e + 2p2 + 0{k\ ky,p^) 

L = [(1 + e)Ji + lOeK]^ + 2[h + + 0{k\ kY^p^) 

Thus, the discontinuity can be measured by 

lim S{0,p) - lim S{k, 0) = - h + 5K- " ^ 



p^o k^o 1 + e 

Of course, this is proportional to K. The dependence on E is captured by the parameter 
e, which also enters into the expressions for the integrals Ji and I2 (see the Appendix for 
a discussion and some characteristic values). Consistent with our expectation, we observe 
that the discontinuity increases monotonically with E, from a limiting value of for 12K < 
(3E ^ 1 to 5.48 for infinite drive. 

Originally obtained from field-theoretic considerations, the form which best displays this 
discontinuity is [p!0|,3. 



sik v) = ne + n^p' 

^ T\\k^ + T^P^ + 0{k\kY,P^) 

near the origin. Here, ny and n± measure the strength of thermal noise in the parallel and 
transverse directions, respectively, while ry and t± are the anisotropic diffusion coefficients. 
Under equilibrium conditions, the FDT enforces the equality 

nil rii 



so that the familiar Ornstein-Zernike form re-emerges. By contrast, we conclude that 
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r± Ty 1 + e 

in the driven case. 

Next, we turn to an estimate for the critical temperature. If we had the exact S{k,p), 
we could identify Tc by its divergence at some point. For the usual system in equilibrium, 
we would look for the divergence of S'(0,0), which is the susceptibility (t^^). However, 
with conserved dynamics at half-filling, this quantity is fixed at zero, so that it is necessary 
to consider limk^o5'(k). In our case, as we just pointed out, this limit depends on the 
direction along which we approach the origin. Thus, we appeal first to the phenomenology, 
i.e., only lim^^o S{0,p) diverges as T — > Tc. In terms of (|12|), only n±/T^ oo. In terms of 
co-operative behavior, this corresponds to an instability against ordering into strips parallel 
to the drive only. 

Now, in a high temperature series, where only a finite number of terms can be computed, 
every partial sum is finite. Instead, the radius of convergence must be estimated. However, 
here we have only one non-trivial term! To make any estimate, we turn to a search for the 
zero of S^^ instead. For the equilibrium case, S~q{k,p) = 1 — 2K{cosk + cosp) + 0{K'^), so 
that this procedure leads to Tc = iJ/ks- Of course, this is the same result, had we expanded 
the exact equation for f3c (e~^^'="^ + 2e~'^^'"^ = 1) in powers of (3cJ and kept only the first 
non-trivial term. Remarkably, this is also the mean field critical temperature. 

Turning to the driven case, we consider an expansion of (given that S = 1): 

S~' = l-^ + OiK') , (13) 

and seek the zero of limp^o S~^{0,p). The result is 

Tc{E) = {J/kn) [5 + h/K] . (14) 

Recall that I2 is proportional to K, so that their ratio is a /3-independent number (Table 1). 

Before discussing the implications of this result, let us identify (up to an overall constant) 
the various parameters in Eqn. (|12D using this approach: 
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r|| = (l + e)(l _ Ji) _ lOKe 

= 2(1 -h- 5K) 
n\\ = 1 + e 
n± = 2 

Referring to the discussion of Ji and I2 given in the Appendix, we see that ry > t± for 
all E > 0. This inequality confirms our choice to identify the critical temperature by the 
vanishing of t±. 

Returning to (p!4D, let us discuss critical temperatures in units of J/ks- First, note that 
I2 is negative and monotonically decreasing in e = exp{—f3E) (cf. Appendix). Therefore, 
Tc{E) increases with E, taking its maximum at infinite E where Tc(oo) = 4.640. For the 
smallest PE, defined through the inequality 12i^ < PE <^ 1, Tc{E) approaches its lowest 
value of 4.0. The agreement with the equilibrium result, Tc(0) = 4, is an artifact of the lowest 



order of the expansion only |13]. Clearly, it is gratifying that even the lowest nontrivial order 
of the high-temperature expansion generates a Tc{E) which increases with E, in qualitative 
agreement with simulation data. Finally, let us consider the quantitative implications of our 
results by focusing on the ratio Tc{E)/Tc{0). For E = 00, our approach yields the value 
1.16, while MC simulations result in 1.40 |jl6|. Thus, the series underestimates this ratio. 



which can be understood as follows. The high temperature series is known to overestimate 



critical temperatures, by underestimating fluctuations. However, it has been argued fll 
that the external drive tends to suppress fluctuations, so that we may expect rc(oo) to be 
less sensitive to the numerical errors introduced by the high temperature series than its 
equilibrium counterpart. In this sense, the series expansion should be "better" for a driven 
system. Indeed, we compare the series result of 4.64 to the simulation value 3.18, finding a 
discrepancy of 46%. In contrast, for the equilibrium case, we have 4.00 and 2.27 respectively, 
showing a much higher discrepancy of 76%. 

Finally, we could hope for better agreement of series and the exact as we move 
into higher spatial dimensions, where fluctuations become less important. The results are 
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certainly encouraging in d = 3. In the series approach, we obtain Tc(0) = 6 and Tc(oo) = 
6.34, signaUing an increase of 6% due to the drive. This is in remarkably good agreement 



with the MC data, which show a 7% increase |T^. Thus, even a low-order calculation can 
produce some quantitatively reliable results. 



IV. CONCLUSIONS 

Within a high temperature series, we have derived the equations for the two-point cor- 
relations of the uniformly driven lattice gas to lowest nontrivial order in j3J, but finite jSE. 
The exact solution of these equations provides a successful qualitative description of two 
central features of our model, namely, the discontinuity singularity of the structure fac- 
tor at the origin, associated with power-law correlations in the disordered phase, and the 
anisotropy in the parallel and transverse diffusion coefficients which controls the onset of 
criticality. Specifically, we observe that the magnitude of the structure factor discontinuity 
increases with E, as a measure of how seriously the FDT is violated in the driven system. 
We demonstrate explicitly that criticality is marked by the vanishing of the transverse dif- 
fusion coefficient, resulting in an estimate for Tf.{E) which increases with E, consistent with 
MC data. On the quantitative side, we argue that fiuctuations, largely neglected in a series 
expansion such as ours, tend to increase the ratio Tc(i?)/Tc(0). In higher spatial dimensions, 
where fiuctuations are less relevant, agreement of series and MC data improves, borne out 
by our results in d = 3. 

Quantitative comparisons aside, the approach presented here provides a convenient an- 
alytic complement to MC simulations, since it gives direct microscopic information about 
effective coarse-grained coupling constants such as and ry which appear in field theories. 
It is computationally simpler than dynamic mean-field theory and easily generalized |13[ to 
higher dimensions or other driven lattice models. Thus, it can help to predict qualitative 
phase diagrams and formulate effective field theories for a wide range of non-equilibrium 
steady states. 
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V. APPENDIX 



Here, we give a few details of our calculations. First, we briefly review the conventions 
of our Fourier transforms and then turn to the evaluation of the integrals contributing to 

Since our lattice consists of discrete points, we will let (x, y) be all pairs of integers. For 
our case, it is most convenient to define the functions Ux,y{k,p): 

Ux,y{k,p) = ^expi{kx + py) 

with continuous k,p E [— vr, vr] . The Ux,y form a complete orthonormal set: 

J2 Ux,y{Kp)Ux,y{k\p') = 5{k - k')5{p - p') 

x,y 

dkdpUx,y{k, p)Ux' ,y'{k, p) = S 

xx' ^yy' 

The Fourier transform is defined in the usual way, 

G{x,y) = J S{k,p) expi{kx + py) 

with inverse 

Sik,p) = ^G{x,y) exp[-i{kx + py)] 

x,y 

Note that S{k,p) is real, since G{x, y) is even in both of its arguments. 
Next, we turn to the integrals contributing to (0). In order to exhibit their properties 
succinctly, it will be helpful to write the anisotropic Laplacian in the form 

A{k,p) = A{1 - cos/c) + B{1 - cosp) 

where the values of interest, A = 2{1 + e) and B = 4 will be inserted at the end. It is then 
easily seen from Eqns. that all integrals are of the general form 

— cos ky{l — cospY 



A{k,p) 
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Clearly, Rij{A,B) = Rji{B,A). Specifically, we need all pairs with i,j = 0,1,2,3 

except (0,0). The calculations are simplified by a series of identities, namely, 



3 
2 



3 
2 



A(fc,p) 
A{k,p) 

(1 — cos k)A{k,p) 

(1 — cosp)A{k,p) 

A(A;,p) 
(1 - cosk)'^A{k,p) 
A(A:,p) 

(1 — cosA;)(l — cosp)A{k,p) 



A{k,p) 
[1 — cosp)'^A{k,p) 



ARio + BRoi 

AR20 + BRii 
ARii + BRq2 
'- AR30 + BR21 

AR21 + BR12 



AR12 + -B-R03 • 



It is thus sufficient to compute Rio only. After performing the elementary integral over p, 
the substitution 1 — cos k = 2t generates, up to prefactors, the integral representation of 



Gauss' hypergeometric function |T^. Defining z = A/B, we obtain 

-Rio 



TxB^^ '2' ' 2' 
^=^° = T5;^^^2'^'2'~'^^ 



We note in passing that all of these integrals can be expressed through elementary functions, 
by reducing the hypergeometric functions down to -F(|, 1; |, — 2;) = 2;"^/^ arctan It is 
now straightforward, if somewhat tedious, to compute l\ and /2, and hence ry and rj_, 
as functions oi z = (1 + e)/2, in the region of interest < e = exp(— < 1. Both 
J's are negative and decrease monotonically. Since the explicit forms are not particularly 
illuminating, we quote a few representative values in Table 1. 
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e = 





0.5 


1 


h/K = 


- 0.844 


-0.931 


-1 


h/K = 


- 0.360 


- 0.761 


-1 



Table 1. 

Characteristic values for Ii and I2, in units of K, for different field strengths. 
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